CREATED USING THE RSC LaTeX PCCP ARTICLE TEMPLATE - SEE www.rsc.org/electronicfiles FOR DETAILS 

ARTICLE TYPE www.rsc.org/xxxxxx | XXXXXXXX 

Procedure to construct a multi-scale coarse-grained model of DNA- 
coated colloids from experimental data 

Bianca M. Mladek*'*'^, Julia Fornleitner'^, Francisco J. Martinez- Veracoechea'', Alexandre Dawid'^, 
and Daan FrenkeF 



o 

(N 



o 



Received Xth XXXXXXXXXX 2010, Accepted Xth XXXXXXXXX 20XX 
First published on the web Xth XXXXXXXXXX 200X 
DOI: 10.1039/bOOOOOOx 



We present a quantitative, multi-scale coarse-grained model of DNA coated colloids. The parameters of this model are transfer- 
able and are solely based on experimental data. As a test case, we focus on nano-sized colloids carrying single- stranded DNA 
strands of length comparable to the colloids' size. We show that in this regime, the common theoretical approach of assuming 
pairwise additivity of the colloidal pair interactions leads to quantitatively and sometimes even qualitatively wrong predictions of 
the phase behaviour of DNA-grafted colloids. Comparing to experimental data, we find that our coarse-grained model correctly 
predicts the equilibrium structure and melting temperature of the formed solids. Due to limited experimental information on the 
persistence length of single- stranded DNA, some quantitative discrepancies are found in the prediction of spatial quantities. With 
the availability of better experimental data, the present approach provides a path for the rational design of DNA-functionalised 
building blocks that can self-assemble in complex, three-dimensional structures. 
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In the pursuit of designing materials that self-assemble into 
specific target structures suitable building blocks have to be 
found with interactions that drive the formation of these struc- 
tures. The availability of such tailor-made nano-structured 
materials could open the way to many interesting applica- 
tions^^. In order to program self-assembly into nano-sized 
building blocks, it is crucial that the interactions between these 
building blocks can be tuned. One class of potential 'pro- 
grammable' building blocks are colloidal particles function- 
alised with polymers. Such particles can be designed in many 
different shapes, ranging in sizes from nm to /im^. More- 
over, the precise choice of their polymeric coating, i.e. type, 
length, flexibility, grafting density and architecture of the 
polymers, allows for additional freedom in tuning the inter- 
actions between the particles. Among such systems, DNA- 
coated colloids (DNACCs) have received special attentionl^El, 
mainly because the technology exists to produce specific DNA 
strands quickly and cheaply. These colloidal particles carry 
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short single- stranded (ss) DNA sequences ("sticky ends") 
connected to inert, grafted polymers ("spacers"). Three- 
dimensional aggregates of such colloids can then be formed 
due to the highly specific and temperature-reversible hybridi- 
sation of complementary sticky ends; these are either carried 
by different species of colloids or are part of so-called linker 
sequences that bridge between different colloids. 

The aggregation behaviour of DNACCs can indeed be in- 
fluenced via the properties of the colloids, their polymeric 
coating, as well as the solution in which the particles are 
immersed. Experiments of nano-^ and micron- sized^^KEl 
DNACCs have shown that self-assembly of simple spatially 
ordered structures, such as bcc or fee crystals, is possible. 
However, applications such as photonic band-gap materials 
would require non-close-packed crystals of low coordination, 
such as the diamond structure and despite recent progress in 
the field the design of arbitrarily complex ordered struc- 
tures is presently still challenging. One crucial factor is that 
DNACCs tend to assemble more readily into amorphous ag- 
gregates than into spatially ordered structures'^. The reason 
is that the attractions between the DNACCs are strongly de- 
pendent on the ex ternal conditions such as ionic strength or 
temperature' ' '. Exquisite control over these parameters is 
thus needed to help DNACCs to anneal into ordered structures. 
Unless we improve our ability to design DNACCs that assem- 
ble readily into the desired target structure, the practical use of 
these building blocks remains limited. 

It is for this reason that the use of coarse-grained models 
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is explored both in theoretical approaches'^ and in com- 
puter simulations'^ These models allow for a fast and 
efficient exploration of new design principles of DNACCs. 
This opens the way to develop strategies for crystals to 
form in broader temperature windows and offer greater 
free dom in the design of DNACCs and the structures they 
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Existing models typic ally range from highly simpli fied on es 
(e.g. lattice modelsl^l^or pair interaction approachesl22S]) 
sophisticated models featuring explicit modelling of the DNA 
hybridisation'^'^ ^'^'^ ^^. In addition, many models exploit 
the elastic properties of the DNA strands: very long strands 
can be described by scaling laws'^'^, while short strands 
of double- stra nded (ds) DNA can be represented as rigid 
rodsEEHmm Many of the existing models are qualitative 
and focus on the generic features of DNACC self-assembly 
- typically, these models do not aim to describe any specific 
DNACC system and hence do not exploit the full available 
experimental information about the building blocks. 

However, for the computer-aided design of DNACCs, quan- 
titative, but computationally tractable models of DNACCs are 
needed. By comparison to the qualitative models mentioned 
above, models of quantitative predictive power are rareEEl^ 
and they are most successful at describing micron-sized col- 
loids covered with short dsDNA strands, where a description 
of the systems via pair interactions determined from simula- 
tions has proven successful. However, the DNACCs that have, 
thus far, shown most promise for crystallisation are the ones 
for which the radius of gyration of the ssDNA strands, Rg, is 
of comparable size to the radius of the colloids. Re- In this 
regime, the modelling strategies that are successful for larger 
colloids cannot be applied: the strands are usually too flexi- 
ble to be approximated as rigid rods but too short for polymer 
scaling laws to apply. 

In a recent Lettei^^ we showed that multi-stage coarse 
graining can be used to describe the phase behaviour of 
nano-sized DNACCs functionalised with ssDNA. The present 
manuscript describes in detail the methodology that we have 
developed to arrive at such a multi-stage coarse-grained 
model. The text is organised as follows: In the first three sec- 
tions, we present three steps of coarse-graining in which we 
identify the key degrees of freedom that determine the phase 
behaviour of DNACCs: we develop our most detailed model 
of DNACCs based on experimental data in Sec. |2] Based on 
simulations of this model, we derive the "core-blob model", 
the second level of coarse-graining, in Sec. [3] This model, in 
turn, allows us to perform the final step of coarse-graining and 
calculate effective interactions (Sec.|4]). The expected reliabil- 
ity of the effective interactions to predict the phase behaviour 
of DNACCs is assessed in Sec. [5] Finally, we calculate the 
phase diagram of the chosen DNACCs within both the core- 
blob model and the effective interaction approach in Sec. [6j 



In the Appendices, we detail technical aspects of the present 
work. 



2 Stage 1: Model with explicit DNA chains 
2.1 General outline 

To develop the most detailed level of description of DNACCs, 
a suitable model for the ssDNA strands tethered to the col- 
loid's surface has to be chosen. ssDNA is not a simple poly- 
mer: it is prone to form hairpin structures and knots. An ac- 
curate description of such substructures could be achieved by 
fully atomistic simulations, which are computationally feasi- 
ble at most for sm all sy stems of DNACCs covered with few, 
short DNA strandsE^Ell. Alternatively, rather detailed, coarse 
grained models of DNA such as developed in Ref . 48" could 
be employed. But while this model makes the detailed study 
of hybridisation between several strands of DNA feasibleS^, 
it is still computationally too time-consuming to be employed 
for a system of hundreds or thousands of colloids, each cov- 
ered with dozens of ssDNAs. Fortunately, the formation of 
ssDNA loops and knots is expected to play a minor role for 
DNACCs where commonly DNAs are chosen that are not self- 
complementary; we thus neglect this effect and model the ss- 
DNA strands as freely jointed, charged chains. This model 
captures the most important contribution to the behaviour of 
the ssDNA strands which stems from the electrostatic repul- 
sion of the DNAs sugar-phosphate backbone In view of 
the high Young's modulus of ssDNA^^, the segments of ev- 
ery freely jointed, charged chain are chosen to have a fixed 
Kuhn length Ixuhn = '^Pss^ where p^s is the persistence length 
of ssDNA. The number of Kuhn segments hk per chain is de- 
termined as 

^contoi 
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Here, /contour = {Nb — l)/?o is the contour length of the ssDNA 
strand, with the number of bases per strand, and bo the in- 
terbase distance in ssDNA. The symbol [. . . J denotes the floor 
function. We stress that since /Kuhn > /^o. each Kuhn segment 
represents several nucleotides. Consequently, our model can- 
not capture the precise base sequence of the ssDNA strands; 
sequence dependent effects, such as base stacking, are cap- 
tured only in an averaged way in the choice of bo (see below). 

The conformation of a freely jointed, charged chain is de- 
fined by the positions {ro , . . . , } of the + 1 vertices of the 
chain. We approximate the continuous charge of the backbone 
by effective charges sitting at each of these vertices; two ver- 
tices / and j at distance rij = \ri — rj\ interact with each other 
via a Debye-Hiickel interaction, <I>^j^, given by 
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where D is the dielectric constant of the solvent and 80 is the 
vacuum permittivity (in SI units). The charge per vertex, is 
approximated as 

^contour 

^=(^^' 

with the effective line charge density of ssDNA, v. Finally, the 
inverse Debye screening length, K, in Eqn.[2]is given by 



(4) 



where Na is Avogadro's number and e is the elementary 
charge. P = l/k^T, where denotes Boltzmann's constant 
and T stands for the temperature (all in SI units). The ionic 
strength / of the solution in which the DNACCs are immersed 
is given in mol/1 and the factor 10^ stems from converting 
mol/1 to SI units. 

The riK segments of each ssDNA strand can be divided 
into two classes: the number of sticky end segments ^se = 
[{Nb^sQ — l)^o/^KuhnJ, calculated from the number of sticky 
end bases A/^,se; and ^sp = — ^se, the number of spacer seg- 
ments. A total of A^str ssDNA strands are then grafted to the 
colloid: we attach the DNA strands at their first vertex ro, 
neglecting in our model the hexane-thiol group by which the 
ssDNA strands are experimentally tethered and which is esti- 
mated to have a end-to-end length of ~ 0.8 nm^^. The col- 
loid, in turn, is modelled as a hard sphere of radius Rc that 
cannot be penetrated by the vertices of the DNA strands (see 
Fig.[T]). In experiments^ colloids are typically maximally 
loaded with DNAs and therefore we assume that the anchoring 
points are uniformly distributed on the surface of the colloid 
and — in accordance with experimental evidence^^^ — that they 
cannot diffuse. 

2.2 Chosen values 

Knowing the experimental conditions under which the ref- 
erence experiments were performed, we can determine the 
values of all variables introduced in the last section. As a 
proof-of-concept of our method, we choose to study system 
V from Ref. 7, where a symmetric, binary mixture of gold 
nano-colloids (labelled A and B) of radius Rc ^ 6 nm were 
studied. All colloids were coated with ~ 60 ssDNA strands 
of Nh = 65 bases, out of which A^^ se = 15 bases constitute the 
sticky end. A and B colloids only differed in their sticky end 
sequences, which were complementary to allow for direct hy- 
bridisation between the unlike colloids. These DNACCs were 
assembled in a solution of 0.01 mol/1 phosphate buffer, 0.2 
mol/1 NaCl at pH = 7.1^. Using the Henderson-Hasselbach 
equationEI, the concentration of all ion species / = 1 , . . . , ^ 
in the solution can be calculated; then, the ionic strength, /, is 

given as, / = ^ £ Cizf, where Zi denotes the charge number of 




Fig. 1 (Colour online) Simulation snapshot of our most detailed 
model of DNACCs: a hard sphere colloid (yellow) of Rc = 6 nm is 
dressed with 60 strands of freely jointed, charged chains of 18 Kuhn 
segments (corresponding to 65 bases; various colours) at the solvent 
conditions given in Ref. 7 and at T = 25 °C The 4 last segments 
model the sticky ends and are plotted in grey. 



ion species /. For the present system^, we find that / = 0.21 
mol/1. 

Next, we need to set the persistence length of ssDNA. While 
the persistence length of dsDNA is well known, its value for 
ssDNA is less well established: a broad range of experimen- 
tally measured values has been published, varying between 
0.75 nm and almost 10 nm^^. Here, we use a value of 0.75 nm, 
in accordance with the studies on which our ssDNA model 
is based According to Ref. \55\ p^s depends on the ionic 
strength and the dependence is approximately described by 
PssiM ^ 4/ [mol/1] Therefore, at the ionic strength of 
our reference experiment (i.e. / = 0.21 mol/1), this results in 
Pss ^ 0.87 nm, which is reasonably close to the value of 0.75 
nm used here. In a similar way, values reported for the in- 
terbase distance bo vary considerably, since they depend on 
the precise DNA sequence under study and the physical con- 
ditions of the solution in which the DNA is immersed. While 
an inter-phosphorus distance of 0.59 nm has been established 
for ssDNA^^, stacking of bases leads to an interbase distance 
that is on average shorteilS^EIl motivated by the findings of 
Ref. [55J we choose a value of bo = 0.43 nm. Therefore, the 
contour length of the ssDNA strands used in the present study 
is /contour ^ 27.5 nm. With both /contour and /xuhn = 2/7ss = 1.5 
nm ready at hand, we find that hk = l^, ^se = 4, and conse- 
quently ^sp = 14. For an isolated DNA strand, we find a radius 
of gyration of Rg ^ 0.1 Rc, thus the radius of gyration of the 
DNA strands is indeed of the order of the colloidal size. 
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Fig. 2 The height distribution //(r, of reactive ends, as function of 

(a) the distance r of the centre of mass of the last 2«se Kuhn segments 
of a strand, Tgnd^ from its anchoring point on the colloid, Fanchor^ and 

(b) the deviation in angle ^ between the vectors Tgnd — ^anchor and the 
connection vector from the colloid's centre to the anchoring point. 
This distribution is peaked at (r ~ l.2Rc, ^ = 0). 

The effective line charge density V is interpolated from Ta- 
ble 1 in Ref. 50 and we find v(/ = 0.21 mol/1; ssDNA) = 2.07 
e/nm, translating to a charge of ^ ~ 3 ^ per vertex. In our study, 
we use Z) = 80 for all temperatures and the Debye screening 
length varies from 0.67 nm (at 25 °C) to 0.71 nm (at 65 °C). 

2.3 Simulations 

To study the behaviour of an isolated DNACC, we implement 
Monte Carlo simulations utilising crankshaft and pivot moves 
to equilibrate the ssDNA chains. Every few Monte Carlo 
sweeps, we also try to regrow whole chains by employing con- 
figurational bias Monte Carlo simulations^^. 

These simulations then allow us to gain insight into the 
height distribution H{r^'6) of the ends of the DNA strands. 
For convenience, we measure this distribution as a function 
of two parameters: (a) the distance r of the centre of mass of 
the last 2^se Kuhn segments of a strand. Tend, from its anchor- 
ing point on the colloid, Tanchor- The choice of 2^se will be 
motivated in Sec. [3j (b) the deviation in angle i^ between the 
vectors Tend — Tanchor and the connection vector from the col- 
loid's centre to the anchoring point. This distribution captures 
how sticky ends are restricted in their movement due to the 
fact that the ssDNA strands are tethered to the colloid and due 
to neighbouring DNA strands, showing a peak at (r ~ 1.27?c, 
1^ = 0) (Fig.[2]). 

In addition, we use a modification of Widom's particle in- 
sertion (mWPI) technique to determine the steric repulsion 
^2,rep(^7 T) between two DNACCs separated by distance r and 
at temperature T in the zero density limit. We find P<I>2,rep to 
be temperature-independent over a wide range of temperatures 
from T=25 °C to 75 °C (Fig. [3] soHd line). 

The large number of degrees of freedom with which 




Fig. 3 The steric repulsion between two DNACCs, |3<l>2,rep^ as func- 
tion of the distance between them and as obtained by the model with 
explicit DNA chains (solid line) and the core-blob model (without 
correction: dash-dotted line, with correction: dashed line). The re- 
pulsion is found to be independent of temperature. The inset shows a 
close-up of the functions for distances r/Rc^5. 

DNACCs are described in the present model would render 
large-scale simulations of crystals of DNACCs unfeasible. 
Also, the present model would call for a binding scheme of the 
sticky ends, w here sev eral Kuhn segments align to form the 
dsDNA stretchE^EHSHSl^ which is computationally rather ex- 
pensive. We therefore refrain from implementing binding be- 
tween complementary ssDNA sequences in the present model 
and rather develop a more coarse-grained model — which we 
term core-blob model — in the next section. 

3 Stage 2: Core-blob model 

The aim of the core-blob model is to arrive at the simplest 
possible model of DNACCs which preserves (i) the steric 
repulsion <I>2,rep between two isolated DNACCs and (ii) the 
height profile of sticky ends with respect to their colloid as 
obtained from the model of explicit DNA chains (see Sec. [2]). 
We therefore model each of the Nstr sticky ends as an entity 
called "blob", which is grafted to the surface of the colloid at 
fixed anchoring points. In this, the last 2^se segments of every 
freely jointed, charged chain constitute a blob; then — to a first 
approximation — the blob's centre represents the connection 
point between the sticky end and its spacer. The gold colloid 
and the remaining hk — 2^se segments of all A^str strands form 
the "core", leading to a model of A^str + 1 separate entities. The 
model is defined by four different interactions (Fig. [4]) which 
we derive from the model of explicit DNA chains via Monte 
Carlo simulations using mWPI^^ and biased simulations'^ 
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(i) the repulsive interaction <I>bb {r) acting between two blobs 
tethered to different colloids and a distance r apart is 
approximated as the interaction of two isolated (=non- 
tethered) freely jointed, charged chains of length 2^se^^- 
As anticipated from studies of polymers, this potential is 
of Gaussian shape (see e.g.^H). However, due to the finite 
length and the charge carried by the ssDNA, the repulsion 
found here is considerably stronger than the 2kBT char- 
acteristic of polymers in the scaling regime (see Fig.|4^); 

(ii) the interaction <l>cb(^) of a blob with the core of another 
colloid separated by distance r. We approximate this 
potential by simulating a single freely jointed, charged 
chain of 2^se segments interacting with a bare core, i.e. a 
hard sphere grafted with chains of length nx — 2^se (see 
Fig#); 

(iii) the interaction ^cc(^) between two cores at distance r 
is estimated as the zero density repulsion between two 
colloids each grafted with chains of length nx — 2^se (see 
Fig.|5); 

(iv) the interaction of a single sticky end with all the re- 
mains of its own colloid (i.e. the core and all other blobs) 
cannot trivially be split into a repulsive and a tether- 
ing contribution due to intricate multi-body contribu- 
tions. In an isolated DNACC, we therefore determine 
the full multi-body potential <I>scb(^5''^) from the height 
distribution of sticky ends in the underlying model as 
^scb(^,''^) = — log//(r,i^) (see Fig.[4]l). By construction, 
this potential guarantees the preservation of the height 
profile of sticky ends with respect to the model of explicit 
DNA chains in the regime of dilute solutions of DNACCs 
(see Sec. [2]). 

To evaluate the reliability of the core-blob model, we de- 
termine the steric repulsion <I>2,rep(^) between two isolated 
DNACCs within this model via mWPI^o and compare the re- 
sults to our findings from the model of explicit DNA chains. 
We find that the core-blob model underestimates <I>2,rep(^) (see 
Fig. [3] dash-dotted line). The reason for this discrepancy can 
be traced back 'to the fact that the core-core repulsion is 
too soft (Fig. [4]:, dashed line), since this potential should also 
include a multi-body contribution from the sticky ends teth- 
ered to the spacer chains, which cannot be captured by <I>bb 
alone. We therefore introduce a correction to <I>cc (Fig. [4]:, 
solid line) by simulating the repulsion between two colloids 
dressed with chains of length nK — ^se instead of nK — 2^se- 
Then, the core-blob model recovers <I>2,rep with sufficient ac- 
curacy for all distances. Especially, we find good agreement 
for distances r > 5.5Rc (see Fig. [3] dashed line) which in- 
cludes the range of experimentally observed next neighbour 
distances, r = aV3/2 > 6.1Rc, where the measured CsCl lat- 
tice spacings a > 42.5 nw^. Defining the colloidal packing 




Fig. 4 (Colour online) The core-blob model, in which sticky ends are 
modelled as blobs and the spacers and colloids constitute a core [see 
inset in a)], a) Blobs on different colloids interact via ^^hh'^ b) a blob 
with the core of another colloid via ^^ch'^ c) different cores via ^^cc 
(dashed line: uncorrected potential; solid line: corrected potential); 
and d) blobs with their own colloid via ^^scb (Isolines shown every 
2.5 ksT). Further, blobs on colloids of different identity can bind and 
form short stretches of dsDNA [sketched as blue rods in the inset in 
a)]. 

fraction r| = ^R^N/V (with V the volume of the conven- 
tional unit cell and N the number of DNACCs in this unit cell) 
and assuming the experimentally observed CsCl structure to 
be the thermodynamically stable structure, we can determine 
the packing fraction up to which our model shows high re- 
liability as Tj = ^Rl2/a^. Using a = 25.5/V3Rc, we find 
11 < 0.033. 

Finally, the hybridisation of complimentary sticky ends has 
to be modelled via a suitable Monte Carlo move. In this, 
we have to account for binding of initially unbound sticky 
ends, breakage of initially bound sticky ends, as well as for 
the change of binding partner for an already hybridised sticky 
end. Further, we wish to use the experimentally measured data 
on the DNA hybridisation free energy. Since the persistence 
length of dsDNA far exceeds that of ssDNA, we model the 
hybridised sticky ends as a (volumeless) rigid rod of fixed 
length L = ^se^Kuhn- For simplicity, we ignore the change 
in inter-base distance between ssDNA (0.43 nm) and dsDNA 
(0.34 nm). Binding is possible between a chosen (bound or 
unbound) blob / on a given colloid and all unbound blobs j 
tethered to unlike colloids and within a distance of approach 
rij < L. Upon binding, the reaction partner j is moved to a 
distance L from / along the connection line rtj. Since bound 
blobs cannot move independently anymore, binding leads to 
the loss of a degree of freedom which is reintroduced upon un- 
binding by placing j along Yij with a probability of r?-, thereby 
guaranteeing detailed balance. The probabilities for each pos- 
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Table 1 The values of the hybridisation free energy Ghyb [^b of the 
DNA strands used in system V in Ref . 7 at different temperatures (ac- 
cording to DINAMelt^^): the DNA spacer sequence is given as 5'- 

AGC TGT GTT CT-3'. Sticky ends on A-colloids read 5'-TAA CCT 
AAC CTT CAT-3', while on 5-colloids the complementary sequence 
is found, 5'-ATG AAG GTT AGG TTA-3'. Values were determined 
at a ionic strength of 0.21 mol/1. 



sible bound state ij and the unbound state / to occur are deter- 
mined by their respective weights W \ 



K 



and 



(5) 



(6) 



where po is the standard density of 1 mol/1, Uij is the potential 
energy of the state where / is bound to j and Uf the potential 
energy of the state where / is unbound; K is the equilibrium 
binding constant of two sticky ends and is connected to the 
hybridisation free energy Ghyb via K = exp (— PGhyb); ^hyb 
(and thereby K) depends on the nucleotide sequence of the 
sticky ends and it is temperature- and salt-dependent; it can be 
approximated as the hybridisation free energy of two sticky 
ends free in solution (e.g. via DINAMelt^^). The values for 
Ghyb used here are given in Tab. [T] The present binding move 
is justified in more detail in the Appendix A. 



4 Stage 3: Effective interactions 

We use the core-blob model to calculate the pair interactions 
between two DNACCs in the zero density limit. Interactions 
between like colloids (i.e. AA and BB) are described by the 
purely repulsive, temperature-independent steric repulsion be- 
tween two DNACCs, P<I>2,rep. Calculated above. The interac- 
tion of unlike colloids (i.e. AB) additionally features an at- 
tractive potential ^2,hyb stemming from hybridisation of DNA 
strands. Thus, 



and 

|3<Df (r, T) = |3<I>2,rep(r) + P<l>2,hyb(r, T). (8) 
P^2,hyb can be obtained by evaluating 

poo 

P*2,hyb(r, T) = - ^ d(|3G(,yb) m)K=.M-PG',^,) ' 

as described in Ref. 32 Here, ^(r) is the number of DNA 
bridges formed between the two colloids at fixed distance r, 
and (...) denotes the statistical average. In practise, the in- 
tegration is performed between Ghyb of interest and a value 
sufficiently large for sticky ends not to hybridise anymore. 

Due to the temperature-dependence of Ghyb, the depth of the 
minimum of ^^{r^T) varies strongly with temperature. For 
instance, a change in temperature from 62.1 to 56.9 °C results 
in a drop in the minimum of <I>^^ of roughly 20 (cor- 
responding to 13 kcal/mol). This strong temperature depen- 
dence of the DNA-mediated attraction explains the difficulty 
in crystallising the DNACCs: only in a narrow temperature 
range around ~ 65°C is the minimum in <$^^{r^T) shallow 
enough to allow for the formation and breakage of DNA links. 
Upon lowering T, the bonds that form cannot break anymore 
and the system gets stuck in disordered aggregates, even if an 
ordered structure is thermodynamically stablel^. 

To assess the predictive power of the core-blob model, it 
would be desirable to compare and to experimen- 

tal results. But while such potentials can be experimentally 
determined for m/cr6>^- sized colloids by using optical tweez- 
gj,gi2 43^ the same is not feasible for nano-colloids. Other ex- 
perimental validation techniques will be required. Thus, val- 
idation of our model against experimental data will only be 
studied later, by comparing the computed and experimentally 
determined phase diagrams. 

Thus far, our approach has allowed us to compute the effec- 
tive pair potential between DNACCs. Pairwise additive inter- 
actions are typically used to model DNACCs as structureless 
particles in theoretical studies (e.g. Ref. ^) and also in some 
computational studies (e.g. Ref. [33]). However, in the regime 



P<^>2,rep(r) 



(7) 



where Rg ~ Rc we may expect that the assumption of pair- 
wise additivity breaks down. With the present model we can 
quantify the importance of such many -body interactions. 



5 Three-body interactions 

Using the core-blob model we test if the three-body interac- 
tions of a system of two A colloids and one B colloid can 
be written as the sum of the various two-body contributions. 
It can be anticipated that differences between the three-body 
interaction and the sum of the two-body contributions will 
mainly arise from a competition of the two A colloids for the 
sticky ends of the B colloid and will therefore crucially de- 
pend on the arrangement of the DNACCs with respect to each 
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Fig. 5 (Colour online) The effective interaction between unlike col- 
loids P^fep (bold lines) is the sum of the temperature-independent 

steric repulsion P<I>2,rep(= P<l>^^^^) between two colloids, and a 
temperature-dependent attractive potential P<J>2,hyb arising from hy- 
bridisation of DNA strands. Data are shown for different tempera- 
tures (top to bottom: 65.1°C, 63.2°C, 6L4°C, 59.1°C, and 56.9°C). 
The inset shows a simulation snapshot of an A colloid (green) inter- 
acting with a B colloid (red). Unhybridised sticky ends are shown 
as yellow spheres, while DNA bridges are shown as blue rods. The 
translucent spheres indicate the average position of the sticky ends. 

other. A linear arrangement of the colloids, with the B colloid 
positioned between the two A colloids, is expected to lead to 
little discrepancy since in the regime studied here (Rg ^ Rc), 
the DNA strands are too short to reach complementary sticky 
ends at the back of the DNACC they face. 

However, in the typical arrangements of DNACCs that oc- 
cur in a crystalline environment, many-body effects are more 
likely to arise. In a crystal, a given colloid is typically sur- 
rounded by several DNACCs of the other species that all com- 
pete for the strands of the central colloid. It is therefore in- 
teresting to study the three-body interaction c^^"^ = [<I>f'^^] 
for three DNACCs arranged in an equilateral triangle of side- 
length r and test if the following relation holds 

p<I>f^^(r,r) = 3P<D2,rep(r) + 2p<I>2,hyb(r,r). (10) 
5.1 Methods 

In analogy to the AB two-body interactions, also the ABA /BAB 
three-body interaction can be split into a contribution <I>3,rep 
stemming from steric repulsions between the DNACCs, and 
an attractive part <I>3,hyb arising from DNA strand hybridisa- 
tion, i.e. 

p<I>f^^(r,r) = p<D3,rep(r,r) + pcI>3,hyb(r,r). (11) 



We first calculate <I>3,rep by generalising the mWPf^ we 
place three non-reactive DNACCs, i.e. DNACCs with non- 
binding sticky ends, in an equilateral triangle of side-length 
^ = ^max bigger than the expected range of the interactions. 
Then we repeatedly reduce the side-length to r — Ar and mea- 
sure (exp [— PAf/ (r ^ r — Ar)]), where MJ {r ^ r — Ar) is the 
change in potential energy between the three DNACCs due to 
the move. In this way, we move the DNACCs together. Then, 
P<I>3,rep(r) = - i:;tt^ log (exp [-PAC/(r' ^ / - Ar)]). 

To determine <I>3,hyb(^). we arrange the three DNACCs in 
an equilateral triangle of fixed side-length r. Measuring the 
total number of DNA bridges formed at this distance, ^(r), we 
can calculate <I^3,hyb(^) as 

poo 

P<I>3,hyb(r) = - J^^^^^ d (PG(,yb) (C(r))^=exp(-K,,) ' ^12) 

in analogy to the determination of ^2,hyb (see Eqn.[9]). 
5.2 Results 

As can be seen from Fig. [6j the repulsive three-body poten- 
tial P^3,rep(^) is, to a good approximation, equal to the sum 
of the two-body contributions 3P<I>2,rep(^)- However, pairwise 
additivity does not hold for the attractive part of the three- 
body potential, P<I>3,hyb- As anticipated, the two A colloids 
increasingly compete for the available sticky ends of B as the 
colloids are moved closer together and fewer bonds can form 
for each of the two AB pairs than in an isolated, single AB 
pair. This overestimation of formed DNA bridges within the 
pair potential approach directly translates to an overestimate 
of the depth of the attraction between the colloids (cf . Eqn. [9] 
and see Fig. [6]) and consequently an underestimate of the po- 
sition of the minimum in the attraction. We therefore expect 
that an analysis of DNACC crystals within the pair potential 
framework will predict more compact crystals than found ex- 
perimentally ^ . 

6 Phase behaviour of DNA coated colloids 

To assess the predictive power of both the core-blob model 
and the pair potential approach, we study the phase behaviour 
of DNACCs by implementing free-energy calculations within 
both approaches. The results are then compared to data avail- 
able from experiments^, such as the stable crystal structure, its 
lattice constant and the melting temperature of these crystals. 

6.1 Crystal structure prediction 

As was already known to OstwaldlSl, observing spontaneous 
formation of a crystal does not imply that the observed struc- 
ture has the lowest free energy. Rather, we have to con- 
sider the thermodynamic stability of all possible crystal struc- 
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Fig. 6 (Colour online) The three-body effective interaction P<l>3(r) 
(bold dashed lines) of two A and one B colloids arranged in an equi- 
lateral triangle compared with the sum of the two-body contributions 
3p<I>2,rep + 2p<l>2,hyb (bold solid lines) at 65.19°C, 63.29°C, 6L49°C, 
59.1°C, and 56.9°C (top to bottom). Further, the repulsive contri- 
bution and the attractive hybridisation contributions to the effective 
interactions are shown explicitly in thinner lines (two-body contribu- 
tions: solid lines, three-body contributions: dashed lines). The inset 
shows a simulation snapshot. A colloids are shown as green spheres, 
and the B colloid as a red sphere. Unhybridised sticky ends are shown 
as yellow spheres, while DNA bridges are sketched as blue rods. The 
translucent spheres indicate the average position of the sticky ends. 



tures. To identify credible candidates for the most stable crys- 
tal structure, w e use optimisation techniques based on ge- 
netic alg orithmsEEH. The structures that the genetic algo- 
rithm identifies as plausible are then considered in the free- 
energy calcul ations . We adapt a search strategy for 2D bi- 
nary mixtureslSSSI to 3D, augmenting it with a parametrisa- 
tion of search space that excludes a priori configurations with 
overlapping colloids . The lattice parameters describing the 
crystal structures are encoded in binary individuals and a ran- 
dom crossover is employed as mating scheme. Mutations take 
place with a rate of 0.05. We limit our search to symmetric 
A5-mixtures and lattice structures with up to eight particles 
per unit cell. Particles interact via their respective pair inter- 
actions (see Sec. [3]), which we first fit to analytical functions 
(see Appendix B). Calculations are run at constant pressure P, 
so that the volume fraction r\ enters the optimisation as an in- 
dependent parameter. Structures are optimised with respect to 
the Gibbs free energy of the system G = U ^ PV — TS, with 
U the internal energy, V the volume of the system, and S the 
entropy. In the genetic algorithm, entropy is neglected, and 
hence the Gibbs free energy is equal to the enthalpy of the 
structures. Temperature-dependence of the system only enters 
our approach via the temperature-dependence of the pair po- 
tential; this limitation in treating entropic effects necessitates 



the subsequent free energy calculations (see Sec. |6.2| ). To de- 
termine the minimum enthalpy configuration at a given pres- 
sure we evaluate 1000 generations of a population of 50 indi- 
vidual ordered structures each. Details on the general working 
principles of the method can be found in Ref . ,65, 

The genetic algorithm calculations predict the CsCl (B2) 
structure as the most stable one for low colloidal packing 
fractions r|, while it also predicts a competing NaTl (B32) 
structure for higher r|. The latter structure is of general in- 
terest since it is composed of two interpenetrating diamond 
structures. If it were possible to remove one of the two col- 
loidal species in a post-assembly modification step, a diamond 
structure could be created Such assembly strategies are ex- 
plored by e.g. substituting the gold colloids in one of the two 
species of DNACCs by organic compounds However, in 
the present system an experimental distinction between CsCl 
and NaTl structures would prove challenging: here, A and B 
colloids only differ by their sticky ends while X-ray scatter- 
ing only detects the gold colloids. As a result, both CsCl and 
NaTl structures would experimentally be detected as bcc ar- 
rangements. 



Apart from the CsCl and NaTl structures, we chose to con- 
sider a few more candidate structures: CuAu (LIq), NaCl 
(Bl), 'straight' hep (s-hcp)^^ ZnS (B3; diamond)"^, AuCd 
(B19), as well as substitutionally disordered CsCl and CuAu 



crystals'^ (Fig. (t)). 
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Fig. 7 (Colour online) The different crystal structures considered in 
this study, a) CsCl; b) NaTl; c) CuAu; d) NaCl; e) ZnS; f) 'straight' 
hep; and g) AuCd. Different colloidal species are coloured red and 
green; for clarity, the DNA strands have been omitted. 



6.2 Free energy simulations 

To study the phase behaviour of DNACCs and to determine the 
stable crystal structure, we have to determine the free energies 
F of all candidate structures mentioned above. Within the ef- 
fective pair potential framework, we calculate F via ther mody- 
namic integration in the canonical ensemble as detailed inl^^^ 
and using systems of at least N = 1000 DNACCs. Within the 
core-blob model, the thermodynamic integration is achieved 
in two steps similar to the approach followed in Ref . [30l using 
systems of more than N = 100 DNACCs. In the first step of 
the integration we arrange N DNACCs in the desired crystal 
structure at fixed colloidal volume fraction r\ and temperature 
T (and thereby fixed binding free energy Ghyb; see Tab. [T]). 
We then transform this DNACC crystal into a system of non- 
binding DNACCs artificially fixed to their lattice sites. We 
achieve this by gradually increasing Ghyb so that sticky ends 
do not bind anymore; at the same time, we gradually confine 
the centres of the colloids to individual, small cells of vol- 
ume V around the perfect lattice sites R^^ of the chosen crystal 
structure by raising a potential barrier. For a colloid / centred 
at Rf and at integration point ^ G [0, 1], the potential energy 
^barr barrier is given as 

[Apt/barr elsC , 

where LJjv(Ry^) is the union of the confining volumes around 
all perfect lattice sites R^^ and (7barr is the maximal height of 
the barrier, t/barr is chosen sufficiently high for the crystals 
not to melt during the thermodynamic integration. The total 
potential energy of the system due to the barrier is given as 
P^barr = ^4=1 Wbarr(^f )' Gauss-Lobatto quadrature^^. 



we numerically evaluate 



|3AFi=|dX(|3<,)- I d(PG0hyb(O^=exp(-K,b)' 

1 PGhyb 

(14) 

where AFi is the difference in free energy between the crys- 
tal of interest and the crystal of inert, confined DNACCs. In 
practise, the upper limit of the second integral is replaced by 
a hybridisation free energy sufficiently large to guarantee that 
no hybridisation of the ssNDA sticky ends takes place. ^ de- 
notes the total number of DNA bridges formed in the system. 

In the second integration step, we use the lattice-coupling 
expansion method we linearly expand both the crystal of 
inert DNACCs and the potential barrier — which guarantees to 
hold the crystal in place — i.e. 

Rf^yRf (15) 

with X = C, LS and the expansion factor y G [1 , . In practise, 
the infinite expansion is approximated by expanding the sys- 
tem sufficiently for particles not to interact anymore. The free 
energy difference between the unexpanded and the expanded 
crystal of inert DNACCs is given by 

PAF2 = p| dy(^)^, (16) 
1 

where "M^ is a modified virial given as 

i<j 

Here, ffj^y is the force acting between DNACC / and DNACC j 
at expansion factor y, and R^- ^ = R^^ — Rf^ is the separation 
of the two DNACCs at the original colloidal packing fraction 
r| (i.e. aty= 1)^^. 

The total free energy per colloid is then given as 

^F/N = P/id + PAFi/A^ + PAF2/A^, (18) 

where P/id is the free energy of an isolated DNACC; being 
the same for all crystal structures, this contribution can be ne- 
glected in the determination of the thermodynamically stable 
crystal structure. 

6.3 Results 

For all candidate crystal structures, we determine the free en- 
ergy for a range of packing fractions, concentrating on the 
regime where T < 65 °C. This is the temperature where the 
effective interaction <I>2 develops a minimum indicating that 
at these temperatures DNA bridges form between DNACCs. 
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Fig. 8 (Colour online) The excess free energy Pr\{F/N — f^^) as 
function of the colloidal volume fraction r| at T = 56.9°C. Both the 
core-blob model (solid lines) and the pair potential approach (dashed 
lines) find the CsCl structure (A) to be the most stable one. This 
crystal is in equilibrium with a dilute vapour (A). The common tan- 
gents are shown as dotted lines. In contrast, the liquid (^) is found 
to be only metastable. In close competition with the CsCl structure, 
a range of metastable crystal structures is found: s-hcp: •, CuAu: ■, 
NaTl: disordered CsCl: *, disordered CuAu: x. Inset: A com- 
parison of the lattice constant a (in nm) as function of temperature 
T (in °C) as measured in the experiments (heating: A, cooling: T^) 
and as obtained from simulations with the core-blob model (•) and 
using the pair potential approach (■). The shaded region indicates 
the temperatures at which crystals are not stable in experiments. 



In this regime, crystals of low packing fraction are expected to 
be stabilised by DNA hybridisation. 

We first present results obtained via the core-blob model 
(Fig. [8] solid lines), initially concentrating on a temperature 
ofr = 56.9°C. We find that at all packing fractions r\ consid- 
ered, the CsCl structure is the most stable, in agreement with 
experimental findings^. In close competition with the CsCl 
structure, we find a whole range of metastable crystal struc- 
tures, namely s-hcp, CuAu and NaTl. These structures are me- 
chanically stable for packing fractions of r| > 0.026 — 0.035; 
by contrast, the CsCl structure is already found stable for 
r| > 0.023. At high r| 0.065, a metastable AuCd phase 
appears. By contrast, both the NaCl and ZnS structures are 
found to be mechanically unstable and melt at all volume frac- 
tions considered. Extrapolation of the collected data suggests 
that the NaTl structure would become the stable structure for 
r| > 0.07. The reason for this transition from CsCl to NaTl can 
be understood from studying the mechanism stabilising the 
latter structure: we found that NaTl crystals are mechanically 
stable only at packing fractions r| > 0.035, where colloids not 
only form DNA bridges with next neighbour colloids, but also 
with the unlike colloids found in the second coordination shell. 
By contrast, CsCl cannot form DNA bridges the the next near- 
est neighbours since they are all like colloids. At sufficently 
high packing fractions beyond r| ~ 0.05, the excluded volume 
effects of the DNA strands render binding to the nearest neigh- 
bours incresingly challenging for both CsCl and NaTl. While 
this severely limits the amount of possible DNA bridges for 
CsCl, NaTl can instead bind to the next nearest neighbours, al- 
lowing it to eventually become the thermodynamically stable 
structure around r\ > 0.07. However, such high volume frac- 
tions cannot be achieved in experiments where crystals form 
from a dilute vapour. Moreover, we stress that the extrapola- 
tion should be taken with a grain of salt as we do not expect 
our model to be fully valid at these high densities (see Sec. [3]). 

Next, we study substitutionally disordered crystals. At low 
volume fractions r\ < 0.05, such crystals have only limited 
mechanical stability and few substitutional changes can be 
sustained by the crystals. A single substitutional defect, in 
which only a single pair of neighbouring A and B colloids 
exchange sites in an otherwise perfect crystal increases the 
free energy of the crystal by AF = OMbT (r| = 0.024) to 
IMbT (r| = 0.05) for CsCl and by AF = OARbT (r| = 0.026) 
to 1 .O^B T(r[= 0.05) for CuAu. At higher r| - 0.05, more sub- 
stitutional disorder can be stabilised: however, the free energy 
of the substitutionally disordered structures is higher than that 
of the perfectly ordered crystals. This can be understood by 
considering the contributions to the free energy. While ^AF2 
(Eqn. 16) and the first integral of ^AFi (Eqn. 14) are inde- 



pendent of the distribution of the A and B colloids over lattice 
sites, the second integral determining ^AFi (Eqn. [14]) depends 
crucially on the average number of DNA bridges {Q formed 
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in the crystal. As {Q is lower in a substitutionally disordered 
crystal than in a perfectly ordered crystal, ^AFi is lower for 
the latter systems. Still, substitutionally disordered crystals 
have been seen to form spontaneously in simulations^^. The 
present results suggest that these structures are kinetically ar- 
rested for r < 65 °C. 

At high r| > 0.1, it is to be anticipated that the phase be- 
haviour of DNACCs will increasingly be dominated by the 
steric repulsion from the compressed DNA strands as well as 
excluded volume interactions of colloids. Then, close-packed 
structures such as CuAu or s-hcp are expected to be thermo- 
dynamically stable. Since crystals are not stabilised by DNA 
hybridisation anymore in this regime, substitutionally disor- 
dered crystals should be favoured over ordered structures. We 
emphasise that exploration of this regime is beyond the scope 
of the current contribution. 

Compared to the results of the core-blob model, the pair 
potential approach — as expected — underestimates the free en- 
ergy by overestimating the number of DNA bridges formed in 
the system (Fig. [8j dashed lines). Still, it offers a good esti- 
mate of the range of mechanical stability of the various crys- 
tal structures and predicts the same phase order as the core- 
blob model, i.e. CsCl as most stable structure, followed by 
metas table s-hcp, CuAu and NaTl. However, within the pair 
potential approach, the NaTl structure is found to out-compete 
the CsCl structure already around lower packing fractions of 
r\ > 0.055. Structures found to be mechanically unstable in 
the core-blob approach (ZnS, CuAu) are also found unstable 
within the current framework. Therefore, the pair potential 
approach offers an excellent tool for assessing the mechanical 
stability of the candidate structures and can be used as a pre- 
selection tool for choosing the crystal structures to be studied 
with the core-blob approach. 

The next question is: what happens when crystals of 
DNACCs melt? Do they form a dilute vapour or a dense liq- 
uid? To answer this question, we determined the melting be- 
haviour of the CsCl structure for both models, and for both we 
find that the crystal coexists with a dilute vapour. A similar 
conclusion was reached in Ref . We find no evidence for a 
phase transition between the dilute solution and a denser liquid 
phase. Using the common-tangent construction (see Fig. [5]), 
we can determine the colloidal volume fraction of the CsCl 
structure at coexistence, r\s- Within the core-blob approach, 
we find r[s = 0.029 at T = 56.9 °C. As expected, the pair po- 
tential approach predicts a more compact equilibrium CsCl 
crystal of r\s = 0.036. Determining r\s within both models 
for several temperatures, we can then determine the predic- 
tion of the lattice constant a of the equilibrium CsCl crystals 
as a function of temperature 



The computed values of a{T) are compared to experimental 
data in the inset of Fig. [8] Concentrating first on the core-blob 
approach, we find that it predicts the thermal expansion coeffi- 
cient of the crystals at least qualitatively correctly. It also cor- 
rectly predicts that at sufficiently low temperatures a levels off 
to a constant value. However, the simulations predict denser 
crystals than experimentally observed. The discrepancy in lat- 
tice constant is as large as ~ 12%. This observation points to 
a problem in the input in our model. One drawback is that 
we neglected the hexane-thiol linker grafting the DNA strands 
to the colloids, which has an end-to-end length of ~ 0.8 nm. 
Taking this linker into account is expected to reduce the dis- 
preancy to the experiments slightly. However, the major weak 
spots are the choices of the 'preferred' experimental values for 
the persistence length /?ss and for the interbase distance bo of 
ssDNA. We chose an average value for both parameters but, 
in reality, both numbers are expected to depend on the pre- 
cise base sequence of the ssDNaE3E1. With more systematic 
experimental data on the sequence dependence of these val- 
ues, we expect that the core-blob model would also allow for 
quantitative predictions of spatial quantities (such as the lattice 
constant). 

We note that the pair potential approach captures neither 
the length nor the temperature behaviour of DNACC crystals 
correctly: it seriously overestimates the thermal expansion co- 
efficient of the crystals. Hence, the pair-potential approach 
cannot be used to describe the thermal properties of DNACC 
crystals. 

For completness, we point out that in the temperature 
regime where DNA strands cannot bind anymore (i.e. well 
above 65 °C), crystallisation can only occur due to excluded 
volume effects for which a high osmotic pressure is needed. 
Since the focus of our work was on crystal formation triggered 
by DNA hybridisation, we did not study this regime. 

The equilibrium melting ('sublimation') temperature of a 
DNACC crystal depends on the concentration of the dilute 
solution and, in contact with an infinitely dilute solution, all 
DNACC crystals will eventually evaporate. However, the 
rate at which this happens depends strongly on temperature. 
The more relevant question is therefore: at what temperature 
does the rate of sublimation of a DNACC crystal become ex- 
perimentally observable? Experimentally, the effective melt- 
ing temperature was determined via ultraviolet- visible spec- 
trophotometry in Ref. 7 and found to be = 62.5 (±0.3) 
°C. To be able to compare simulational data to experiments, 
we can estimate the temperature below which the spontaneous 
evaporation of DNACCs from a crystal becomes negligible. 
A rough estimate of the concentration where this happens can 
be obtained using Smoluchowski's treatment of the diffusion 
limited growth of a cluster*^^^. In equilibrium, the evapora- 
tion rate of DNACCs from a solid, spherical cluster of radius 
R equals the aggregation rate to that cluster in the presence of 
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a dilute vapour of DNACCs of density py 
the latter rate, &N /&t, as 



dA^ 
dT 



We can determine 



(20) 



with D being the diffusion constant of a DNACCpSSl. On the 
other hand, the amount of DNACCs in the solid cluster can be 
expressed slsN = ^R^ps, where p^ is the density of the cluster 
and p^ = r\s/{^Rc)- Then, the following relation between py 
and p^ can be derived: 



1 dR^ 
26 "dT 



(21) 



To estimate py, we assume that a solid cluster grows to 1 jum 
in one day (~ lO^s). Further, we estimate the diffusion con- 
stant of a DNACC from the Einstein-Smoluchowski relation, 
D = ksT / {6nf\Rc) . With fj ~ IcP being the viscosity of wa- 
ter, we find that D ^ 4x 10"^cm^/s. Then, py/p^ 10"^. 
From the knowledge of the coexistence densities of the vapour 
and solid phases at several temperatures, we find that the last 
relation is fulfilled for = 63.5 (±0.2) °C in the pair po- 
tential approach and for Tm = 61.9(±0.2) °C in the core-blob 
model (see Fig. [9]), which is in good agreement with the exper- 
imental finding of Tm = 62.5(±0.3) °C. We note that in Ref. 
[37I we estimated the melting temperature within the core-blob 
approach by considering the point where small colloidal crys- 
tals melted on the time scale of a simulation. That approach 
is likely to lead to a higher estimate of the melting temper- 
ature and, indeed, we found that simulated crystals melt for 
T > 64.3(±0.5) °C (see Fig.[9]). 

7 Conclusion 

In this paper we described a method to construct a quantita- 
tive coarse-grained model of DNACCs and compared its re- 
sults to experimental data. For comparison, we also studied 
the predicitve power of a simpler model based on the effective 
interactions between pairs of DNACCs. 

We found that the pair potential approach can be used as a 
qualitative tool, allowing to roughly delimit the range of me- 
chanical stability of DNACCs, and give a coarse estimate of 
the melting temperature. Further, it allows for qualitative in- 
sight into the compression behaviour of DNACC crystals upon 
temperature reduction. This model has the advantage of being 
computationally inexpensive, therefore allowing for fast and 
large-scale testing of DNACC designs. 

However, we showed that for quantitative insight into 
DNACC systems in the regime where the radius of gyration 
of the tethered DNA strands is of the order of the size of the 
(nano-)colloid, an explicit description of sticky ends is needed 
to capture the competition of DNACCs for DNA bridges. We 
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Fig. 9 Top: The ratio between the coexistence densities of the vapour 
and the solid, Pv/p^, as function of the temperature according to sim- 
ulations within the core-blob model (•) and the pair interaction ap- 
proach (■). The solid and dashed lines serve as guides to the eye. The 
threshold value of 10~^ (dotted line) indicates the temperature below 
which spontaneous evaporation of DNACCs from a crystal becomes 
negligible. Middle and bottom: Simulation snapshots of a DNACC 
crystal (core -blob model) at T = 65.1 °C which melted during the 
course of the simulation and of a crystal dXT = 60.2 °C, where the 
evaporation of DNACCs is negligible. A colloids are shown as green 
spheres, B colloids in red. For clarity, only hybridised DNA chains 
are shown. 
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therefore developed a more detailed model, which we termed 
"core-blob model" and which is solely based on experimental 
input. Results from this approach show good agreement with 
experimental data in temperature-dependent quantities. While 
the lattice constant is not captured quantitatively in absolute 
values, the thermal expansion coefficient of crystals is well 
described. We speculate that, once more systematic experi- 
mental data on the persistence length and inter-base distance 
of ssDNA become available, the core-blob model should also 
account for the experimentally observed lattice spacings. We 
could of course have adjusted the model parameters to account 
for the observed lattice spacing, but this would have defeated 
the purpose of the present work, which was to construct a 
model based exclusively available ssDNA data - our model 
contains no a posteriori fitted parameters. 

A potential drawback of the core-blob model is that the in- 
put parameters depend on the ssDNA length and sequence. 
For a given set of input parameters, new simulations are 
needed to redevelop the interaction potentials of the core-blob 
model and to predict the phase behaviour. In some cases, the 
situation may be better. For instance, in e.g. Ref. 8 , differ- 
ent systems were generated by supplementing one choice of 
DNACCs with different linkers. Then, the length of spac- 
ers, their binding strength and the number of reactive ends can 
be tuned in a straightforward way via these linker sequences 
alone. Similarly, our model can be easily generalised to incor- 
porate linkers while reusing the present representation of the 
DNACCs themselves. Furthermore, the core-blob model can 
be adapted to study e.g. systems of more complex coatings, 
asymmetric mixtures, or polydisperse systems, while allowing 
for direct mapping to the corresponding experimental system. 

In summary, a rough scanning of the phase behaviour of 
DNACC designs via the pair potential approach can be used to 
preselect promising DNACC designs, which can then be quan- 
titatively studied in more detailed, but also computationally 
more expensive calculations using the core-blob model. In this 
way, the approach presented here offers a path to computer- 
aided design of suitable DNA-grafted building blocks, advanc- 
ing the efforts of constructing truly complex self-assembling 
structures. 
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9 Appendix A: The binding move in the core- 
blob model 

In the following, we derive the Monte Carlo algorithm of the 
binding move used for hybridisation of sticky ends in the core- 



blob model. In Sec. |9.1| we study the requirement of detailed 
balance for the case of an unbound sticky end binding to a 
complementary, unbound sticky end within reach, i.e. within 
distance of hybridised sticky ends, L. Without loss of gener- 
ality, we name the chosen sticky end a and assume that it is 
attached to a colloid of kind A. Consequently, we term its pos- 
sible binding partner b, which is fixed on a colloid of kind B. 
We then generalise to an arbitrary number of possible binding 



partners (incl. the possibility of a partner change) in Sec. 9.2 



9.1 Detailed balance for one possible binding partner 

To fulfil detailed balance, we need to justify that the flow ^ 
from the configuration where a is unbound ("a") to the state 
where a and h are hybridised ("ab'') is the same as the reverse 
flow, i.e. 

X{a ah) = %{ab a) . (22) 

We can write 

%{a ah) = PaPgQn{ci ab)Ps^cc{^ cib), (23) 

where Pa is the probability of a being unbound, /gen(<2 ^ 
is the probability that the Monte Carlo move hybridises a with 
b and P^cc{a ab) is the probability of accepting this move. 
An analogous formula can be written for %{ab a). 

Next, we write the various terms in Eqn.|23]in terms of the 
coordinates of the chosen sticky end a that we try to bind, r^, 
and the distance of a from its possible binding partner b, rat = 
|r^ — r^|. All other coordinates will be denoted by {rrest}- 

The probability to be in the unbound state is then given by 



-^^-drarl^drabdQ.{drr^si} 



(24) 



where r^^dratdO. is an infinitesimal volume element around 
the location of b. The potential energy of the state where a 
is unbound is given by Ua, while qf^ and q^^^ are the internal 
partition functions of a and b. A^ and A^ denote the respective 
de Broglie wavelengths. 
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When the sticky ends a and b hybridise, we place along 
the connection line Yab at a distance L from r^. The co- 
ordinate of a remains unchanged. Since there is only one 
way to implement this move, the generation probability is 

/gen(<^ Clb) = 1. 

The probability to be in the hybridised state, Pat, where a 
and b are connected by a rod of length X, is given by 



1111 



(25) 



where Uab is the potential energy of the state where a and b are 
hybridised. The internal partition function of the hybridised 
sticky ends which are restrained in their rotational freedom is 
denoted by q^^^'^^^^- Since q^^^'^^^^ is independent of the precise 
orientation of the hybridised sticky ends in space, it can be 
related to the internal partition of a rotationally unrestricted 



dsDNA segment, q^^L via q 



q^^u/^Ti and thus 



Pab = e-^^^^drada{dr,,,t} 



€1 



(26) 



The factor An results from the fact that the connection rod 
between a and b is restrained from its rotational freedom. Fur- 
ther, Aab is the de Broglie wavelength. 

To generate the unbound state from the hybridised one, we 
have to generate a new coordinate for b along the connection 
line of the bound b and a, Yab- Generating the new position 
with a probability proportional to the distance squared, we get 



(27) 



By imposing Eqn. 22 we get the condition for the acceptance 
probabilities for the Monte Carlo move: 



^acc(<^ Clb) 



K 



(28) 



P^,,{ab^a) f X^po 
with MJ = Uab — Ua- In this last equation, we have used thaP 
€b^l^l _ K 



^f<Kb po' 



(29) 



with K being the equilibrium binding constant and po the stan- 
dard density of 1 mol/1. 

9.2 Algorithm for the binding move in the case of several 
possible binding partners 

Here, we outline the algorithm for the Monte Carlo bind- 
ing/partner change move for several possible binding partners. 



1 . Choose a blob (bound or unbound) at random. Without 
loss of generality, we assume that this blob is on a colloid 
of kind A. Therefore, we denote this blob as a and its 
coordinates as r^. 

2a. Find all 7 = 1,...,A/^ unbound blobs on unlike col- 
loids (here: kind B), {bj}, that are within distance X, 
i.e. \rabH I = —^a\ ^ ^ (see Fig. 10 1). Here, r^w is the 
coordinate of the unhybridised possTBle binding partner 
j- 

2b. In case a is initially bound, add its actual binding partner 
to this list; the coordinate of this binding partner is de- 
noted as r^h since it is hybridised with a. We then have a 

total of A^^ + 1 possible binding partners. 

3 a. For all unbound binding partners bj at their unbound 
positions r^^^, do the following: along the connection 
line r^^M , generate the new position r^h as if bj were hy- 
bridised with a (see Fig. [ToJd), i.e. place it at distance L 
from a: 

(30) 



--^^abyl^ab^jl 



3b. If a was hybridised initially, further compute the follow- 
ing for its actual binding partner bj: Along the connec- 
tion line r^^h , randomly generate a new, unbound position 
of bj, Yuu as 

= r« + {vj^h -Y^. (31) 

with random number x G [0, 1). 

After this step, both r^^h and r^^w are known for all possible 
binding partners. 

4a. Calculate the weight of the state where a is not hybridised 
as 

W« = exp(-Pf/^), (32) 

where 

Ua = Y,Uby (33) 

j 

In the last equation, Ub^. is the energy that each of the j 
possible binding partners at their unbound positions have 
with the rest of the system, including the repulsion of 
a. Care has to be taken not to double count interactions 
between the various bj. 

4b. Calculate the weights for each of the possible hybridised 
states. The weight of the state where bj is bound to a is 
given as 

^ (34) 
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Fig. 10 (Colour online) The hybridisation move, a) From a cho- 
sen blob a on colloid A, all possible (unbound) binding partners bj 
are determined which are within a sphere of radius of the hybridised 
sticky ends, L. b) A binding partner, bj, is chosen according to its 
weight (see text). This sticky end is moved from its unbound coordi- 
nates, b^j, to its coordinates in the hybridised state, b^j. c) The DNA 
bridge formed by the hybridised sticky ends is modelled as a rigid 
rod of fixed length (blue rod). Note that the figure disregards the case 
where a is originally bound. 



where 



Uabj = U 



(35) 



where Uiju are the energies that each of the j) possi- 

ble binding partners at their unbound positions have with 
the rest of the system, including the repulsion of a and hj 
at their bound positions. Further, U^y^ is the interaction 

of bj at its bound position with the rest of the system, in- 
cluding a. Again, care has to be taken not to double count 
interactions. 



5. Calculate the sum S over all weights 

j 



(36) 



6. Randomly choose a state according to the probabilities of 
the unbound state 

W 



and the various hybridised states 

_ Watj 
Pabj - ^ ■ 



(38) 



If a hybridised state abj is chosen, a and bj are connected 
by a rod and cannot move independently anymore (see 
Fig.flOb). 
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Appendix B: Fitting functions for effective 
interactions 



To use effective interactions in computer simulations, it is cru- 
cial to find reliable fits for both <l>2,rep and <l>2,hyb- The former 



can be fitted by a sigmoidal curve for large distances r > 5.5Rc 

(39) 



1 +exp 



with fitting parameters ai, Z?i and ci. For distances r < 5.5/?c, 
the potential can be approximated by an exponential function 



P^2!rep(^<5.57?c)=^2exp| 



r 



(40) 



with fitting parameters ^2, b2 and C2. The fits have to be per- 
formed under the boundary condition that the two functions 
need to join smoothly at r = 5.5Rc. 

The attractive potential <l>2,hyb is fitted in two steps: for dis- 
tances where DNA strands can hybridise, i.e. r < 6Rc, the 
potential is approximately linear, <l>2,hyb ^ kx-\-d. For dis- 
tances r > 7.5Rc, ^2,hyb = 0- Then, <l>2,hyb can be fit for each 
temperature of interest as an interpolation between these two 
trends i.e. 



P^2 hyb(^' ^ fi^^d) = ^-^ + ^^^^111 



1 +exp 



a — x 



(41) 



with a = —d/k and fitting parameter x. 

References 

1 G. M. Whitesides and M. Boncheva, Proceedings of the 
National Academy of Sciences, 2002, 99, 4769-4774. 

2 S. C. Glotzer and M. J. Solomon, Nature Mater, 2007, 6, 
557. 

3 N. Geerts and E. Eiser, Soft Matter, 2010, 6, 4647-4660. 

4 C. A. Mirkin, R. L. Letsinger, R. C. Mucic and J. J. 
Storhoff, Nature, 1996, 382, 607. 

5 A. R Alivisatos, K. R Johnsson, X. Peng, T. E. Wilson, 
C. J. Loweth, M. R Bruchez and R G. Schultz, Nature, 
1996, 382, 609. 

6 S. Y. Park, A. K. R. Lytton-Jean, B. Lee, S. Weigand, G. C. 
Schatz and C. A. Mirkin, Nature, 2008, 451, 553. 

7 D. Nykypanchuk, M. M. Maye, D. van der Lelie and 
O. Gang, Nature, 2008, 451, 549. 

8 H. Xiong, D. van der Lelie and O. Gang, Phys. Rev. Lett., 
2009, 102, 015504. 

9 R. J. Macfarlane, M. R. Jones, A. J. Senesi, K. L. Young, 
B. Lee, J. Wu and C. A. Mirkin, Angew. Chem. Int. Ed., 
2010, 49, 4589. 

10 M. M. Maye, M. T. Kumara, D. Nykypanchuk, W. B. Sher- 
man and O. Gang, Nat. NanotechnoL, 2010, 5, 116. 

11 D. Sun and O. Gang, Journal of the American Chemical 
Society, 2011, 133, 5252-5254. 

12 P. L. Biancaniello, A. J. Kim and J. C. Crocker, Phys. Rev. 
Lett., 2005, 94, 058302. 



This journal is ©The Royal Society of Chemistry [year] 



Journal Name, 201 0, [vol], 1 -[Tt] | 1 5 



13 A. J. Kim, P. L. Biancaniello and J. C. Crocker, Langmuir, 
2006, 22, 1991-2001. 

14 M. T. Casey, R. T. Scarlett, W. Benjamin Rogers, I. Jenk- 
ins, T. Sinno and J. C. Crocker, Nat. Commun., 2012, 3, 
1209. 

15 K. M. Ho, C. T. Chan and C. M. Soukoulis, Phys. Rev. 
Lett, 1990, 65,3152-3155. 

16 R Cigler, A. K. R. Lytton-Jean, D. G. Anderson, M. G. 
Finn and S. Y. Park, Nat. Mater., 2010, 9, 918. 

17 R. J. Macfarlane, B. Lee, M. R. Jones, N. Harris, G. C. 
Schatz and C. A. Mirkin, Science, 2011, 334, 204-208. 

18 Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, 
L. Feng, A. D. Hollingsworth, M. Week and D. J. Pine, 
Nature, 2012, 491,51. 

19 F. Varrato, L. Di Michele, M. Belushkin, N. Dorsaz, S. H. 
Nathan, E. Eiser and G. Foffi, Proceedings of the National 
Academy of Sciences, 2012, X, X. 

20 R. Jin, G. Wu, Z. Li, C. A. Mirkin and G. C. Schatz, Jour- 
nal of the American Chemical Society, 2003, 125, 1643. 

21 R. Dreyfus, M. E. Leunissen, R. Sha, A. V. Tkachenko, 
N. C. Seeman, D. J. Pine and P. M. Chaikin, Phys. Rev. E, 
2010, 81, 041404. 

22 A. V. Tkachenko, Phys. Rev Lett., 2002, 89, 148303. 

23 N. A. Licata and A. V. Tkachenko, Phys. Rev. E, 2006, 74, 
041408. 

24 A. V. Tkachenko, Phys. Rev Lett., 2011, 106, 255501. 

25 P. Varilly, S. Angioletti-Uberti, B. M. Mognetti and 
D. Frenkel, The Journal of Chemical Physics, 2012, 137, 
094108. 

26 F. W. Starr and F. Sciortino, Journal of Physics: Con- 
densed Matter, 2006, 18, L347. 

27 B. Bozorgui and D. Frenkel, Phys. Rev. Lett., 2008, 101, 
045701. 

28 W. Dai, C. W. Hsu, F. Sciortino and F. W. Starr, Langmuir, 
2010, 26, 3601-3608. 

29 F. J. Martinez- Veracoechea, B. Bozorgui and D. Frenkel, 
Soft Matter, 2010, 6, 6136-6145. 

30 F. J. Martinez- Veracoechea, B. M. Mladek, A. Tkachenko 
and D. Frenkel, Phys. Rev. Lett, 2011, 107, 045902. 

3 1 R Vargas Lara and R W. Starr, Soft Matter, 201 1 , 7, 2085- 
2093. 

32 M. E. Leunissen and D. Frenkel, /. Chem. Phys., 2011, 
134, 084702. 

33 R. T. Scarlett, M. T. Ung, J. C. Crocker and T. Sinno, Soft 
Matter, 2011,7, 1912-1925. 

34 C. Knorowski, S. Burleigh and A. Travesset, Phys. Rev. 
Lett., 2011, 106,215501. 

35 C. Chi, R Vargas-Lara, A. V. Tkachenko, R W. Starr and 
O. Gding, ACS Nano, 2012, 6, 6793-6802. 

36 T. I. Li, R. Sknepnek, R. J. Macfarlane, C. A. Mirkin and 



M. Olvera de la Cruz, Nano Letters, 2012, 12, 2509-2514. 

37 B. M. Mladek, J. Fomleitner, F. J. Martinez- Veracoechea, 
A. Dawid and D. Frenkel, Phys. Rev. Lett., 2012, 108, 
268301. 

38 B. M. Mognetti, M. E. Leunissen and D. Frenkel, Soft Mat- 
ter, 2012, 8, 2213. 

39 N. A. Licata and A. V. Tkachenko, Phys. Rev. E, 2006, 74, 
040401. 

40 S. H. Tindemans and B. M. Mulder, Phys. Rev E, 2010, 
82, 021404. 

41 S. Angioletti-Uberti, B. M. Mognetti and D. Frenkel, Nat. 
Mater, 2012, 11,518. 

42 M. E. Leunissen, R. Dreyfus, R. Sha, N. C. Seeman and 
P. M. Chaikin, Journal of the American Chemical Society, 
2010, 132, 1903-1913. 

43 W. B. Rogers and J. C. Crocker, Proceedings of the Na- 
tional Academy of Sciences, 2011, 108, 15687. 

44 B. M. Mognetti, P. Varilly, S. Angioletti-Uberti, F. J. 
Martinez- Veracoechea, J. Dobnikar, M. E. Leunissen and 
D. Frenkel, Proceedings of the National Academy of Sci- 
ences, 2012, 109, E378-E379. 

45 W. B. Rogers and J. C. Crocker, Proceedings of the Na- 
tional Academy of Sciences, 2012, 109, E380. 

46 O.-S. Lee and G. C. Schatz, The Journal of Physical 
Chemistry C, 2009, 113, 2316-2321. 

47 V. A. Ngo, R. K. KaHa, A. Nakano and R Vashishta, /. 
Phys. Chem. C, 2012, 116, 19579. 

48 T. E. Ouldridge, A. A. Louis and J. R K. Doye, The Jour- 
nal of Chemical Physics, 2011, 134, 085101. 

49 T. E. Ouldridge, A. A. Louis and J. P. K. Doye, Phys. Rev. 
Lett., 2010, 104, 178101. 

50 Y. Zhang, H. Zhou and Z.-C. Ou-Yang, Biophys. J, 2001, 
81, 1133. 

51 S. J. Hurst, A. K. R. Lytton-Jean and C. A. Mirkin, Ana- 
lytical Chemistry, 2006, 78, 8313-8318. 

52 D. Nykypanchuk, private communication. 

53 L. J. Henderson, Am. J Physiol, 1908, 21, 173. 

54 S. V. Kuznetsov, Y. Shen, A. S. Benight and A. Ansari, 
Biophys. J, 2001, SI, 2^64. 

55 B. Tinland, A. Pluen, J. Sturm and G. Weill, Macro- 
molecules, 1997, 30, 5763. 

56 S. B. Smith, Y. Cui and C. Bustamante, Science, 1996, 
271, 795. 

57 B. M. Mognetti et al, to be pubHshed. 

58 J. B. Mills, E. Vacano and P. J. Hagerman, Journal of 
Molecular Biology, 1999, 285, 245 - 257. 

59 D. Frenkel and B. Smit, Understanding Molecular Simu- 
lation, Academic Press, London, 2nd edn, 2002. 

60 B. M. Mladek and D. Frenkel, Soft Matter, 2011, 7, 1450. 

61 P. G. Bolhuis, A. A. Louis, J.-P. Hansen and E. J. Meijer, 



1 6 I Journal Name, 20 1 0, [vol]. 1 -[TT] 



This journal is ©The Royal Society of Chemistry [year] 



The Journal of Chemical Physics, 2001, 114, 4296. 

62 N. R. Markham and M. Zuker, Nucleic Acids Res., 2005, 
33, W577. 

63 R. A. Van Santen, /. Phys. Chem., 1984, 88, 5768. 

64 J. H. Holland, Adaptation in Natural and Artificial System, 
The University of Michigan Press: Ann Arbor, 1975. 

65 D. Gottwald, G. Kahl and C. N. Likos, /. Chem. Phys., 
2005, 122, 204503. 

66 J. Fornleitner, F. Lo Verso, G. Kahl and C. N. Likos, Soft 
Matter, 2008, 4, 480. 

67 J. Fornleitner, F. Lo Verso, G. Kahl and C. N. Likos, Lang- 
muir, 2009, 25, 7836. 

68 G. J. Pauschenwein, /. Phys. A: Math. Theor, 2009, 42, 
355204. 

69 B. M. Mladek, P. Charbonneau and D. Frenkel, Phys. Rev. 
Lett., 2007, 99, 235702. 

70 B. M. Mladek, R Charbonneau, C. N. Likos, D. Frenkel 
and G. Kahl, Journal of Physics: Condensed Matter, 2008, 
20, 494245. 

7 1 M. Abramowitz and L A. Stegun, Handbook of Mathemat- 
ical Functions with Formulas, Graphs, and Mathematical 
Tables, Dover, New York, ninth Dover printing, tenth GPO 
printing edn, 1964. 

72 E. J. Meijer, D. Frenkel, R. A. LeSar and A. J. C. Ladd, /. 
Chem. Phys., 1990, 92, 7570. 

73 M. von Smoluchowski, Z. Physik. Chem., 1917, 92, 129- 
168. 

74 Diffusion-Limited Reactions, ed. C. H. Bamford, C. Tipper 
and R. G. Compton, Elsevier, 1985, vol. 25, pp. 3 - 46. 

75 B. Bozorgui, PhD thesis. University of Amsterdam, 2008. 

76 C. C. Kerr, S. J. van Albada, C. J. Rennie and R A. Robin- 
son, Clinical Neurophysiology, 2010, 121, 962 - 976. 



This journal is ©The Royal Society of Chemistry [year] 



Journal Name, 201 0, [vol], 1 -[Tt] | 1 7 



